// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/vtk.hh>
#include "parfvdatahandle.hh"
using namespace Dune;

//{evolve_signature_begin}
template<class GridView,class Mapper>
void paraevolve(const GridView& gridView,
            const Mapper& mapper,
            std::vector<double>& c,                
            double dt,
            const std::function<FieldVector<double,GridView::dimension>
            (FieldVector<double,GridView::dimension>)> v,
            const std::function<double (FieldVector<double,GridView::dimension>)> inflow)
//{evolve_signature_end}
{
  //{evolve_init_begin}
  //dim is the dimension of the reference element
  constexpr int dim = GridView::dimension;
  //dimesion of its world space
  constexpr int dimworld = GridView::dimensionworld;
  //the number type used for coordinates
  typedef typename GridView::ctype ct;
  std::vector<double> update(c.size());
  std::fill(update.begin(),update.end(),0.0);

  // check data partitioning
  assert(gridView.overlapSize(0)>0 || (gridView.ghostSize(0)>0));
  //{evolve_init_end}


  //{element_loop_begin}
  for(const auto& element : elements(gridView,Partitions::all))
    {
      // Element geometry
      auto geometry = element.geometry();
      double elementVolume = geometry.volume();
      typename Mapper::Index i = mapper.index(element);
       //{intersection_loop_begin}
       for(const auto& intersection : intersections(gridView,element))
         {
           auto intersectionGeometry = intersection.geometry();
           // Center of intersection in global coordinates
           FieldVector<double,dim> intersectionCenter=intersectionGeometry.center();
           // Center of intersection in local coordinates
           const auto& intersectionReferenceElement = ReferenceElements<double,dim-1>::general(intersection.type());
           FieldVector<double,dim-1> intersectionLocalCenter = intersectionReferenceElement.position(0,0);
           // Velocity at intersection center v_ij
           FieldVector<double,dim> velocity = v(intersectionCenter);
           // Normal vector scaled with intersection area: n_ij |\gamma_ij|
           FieldVector<double,dim> integrationOuterNormal = intersection.integrationOuterNormal(intersectionLocalCenter);
           // compute factor ocurring in flux formula <v_ij, n_ij> |\gamma_ij|
           double intersectionFlow = velocity*integrationOuterNormal;
           // outflow contributions
           update[i]-=c[i]*std::max(0.0,intersectionFlow)/elementVolume;
          // inflow contributions
          if(intersectionFlow<=0)
            {
              // interior intersection
              if(intersection.neighbor())
                {
                  auto j = mapper.index(intersection.outside());
                  update[i] -=c[j]*intersectionFlow/elementVolume;
                }
              //boundary intersection
              if (intersection.boundary())
                update[i]-=inflow(intersectionCenter)*intersectionFlow/elementVolume;
            }
         }
       //{intersection_loop_end}
     }
   //{element_loop_end}

   VectorExchange<Mapper,std::vector<double>> dh(mapper,update);
   gridView.template communicate<VectorExchange<Mapper,std::vector<double>>>(dh,InteriorBorder_All_Interface,ForwardCommunication);
   // auto consistentUpdate = update;
   // VertexDataUpdate<GridView,Vector> matrixDataHandle(gridView,update,consistentUpdate);
   // gridView.communicate(matrixDataHandle,InteriorBorder_InteriorBorder_Interface,ForwardCommunication);


  for(std::size_t i=0; i<c.size(); ++i)
     c[i]+=dt*update[i];  
}

/**
 * This file come from grid-howto and also chapter 3.4 of the dune book.
 **/
int main(int argc, char** argv)
{
  try{
    // Maybe initialize MPI
    MPIHelper& helper = MPIHelper::instance(argc, argv);
    std::cout << "Hello World! This is getting-started-transport-fv." << std::endl;
    if(MPIHelper::isFake)
      std::cout<< "This is a sequential program." << std::endl;
    else
      std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
        <<" processes!"<<std::endl;
    // return 0;
  }
  catch (Exception &e){
    std::cerr << "Dune reported error: " << e << std::endl;
  }
  catch (...){
    std::cerr << "Unknown exception thrown!" << std::endl;
  }

  // {create_grid_being}
  int overlapSize = 1;
  // --- 1d ---  
  // constexpr int dim = 1;
  // using Grid = YaspGrid<dim>;
  // std::bitset<dim> periodic("0");
  // Grid grid({1.0},{80},periodic,overlapSize,MPI_COMM_WORLD);  
  // --- 2d ---
  // constexpr int dim = 2;
  // using Grid = YaspGrid<dim>;
  // std::bitset<dim> periodic("00");
  // Grid grid({1.0,1.0},{80,80},periodic,overlapSize,MPI_COMM_WORLD);  
  // --- 3d ---  
  constexpr int dim = 3;
  using Grid = YaspGrid<dim>;
  std::bitset<dim> periodic("000");
  Grid grid({1.0,1.0,1.0},{40,40,40},periodic,overlapSize,MPI_COMM_WORLD);  

  using GridView = Grid::LeafGridView;
  GridView gridView = grid.leafGridView();
  // {create_grid_end}

  // {create_concentration_begin}
  // Assigns a unique number to each element
  MultipleCodimMultipleGeomTypeMapper<GridView> mapper(gridView, mcmgElementLayout());
  // {create_concentration_end}  

  // Allocate a vetor for the concentration
  std::vector<double> c(mapper.size());
  // {lambda_initial_concentration_begin}
  // norm: one_norm, two_norm
  auto c0 = [] (const FieldVector<double,dim>& x){
              return (x.two_norm()>0.125 && x.two_norm()<0.5) ? 1.0 : 0.0;
            };
  // {lambda_initial_concentration_end}
  // {sample_initial_concentration_begin}
  // Iterate over grid elements and evaluate c0 at element centers
  for (const auto& element : elements(gridView))
    {
      // Get element geometry
      auto geometry = element.geometry();
      auto global = geometry.center();
      c[mapper.index(element)] = c0(global);
    }
  // {sample_initial_concentration_end}

  // Construct VTK writer
  // {construct_vtk_writer_begin}
  auto vtkWriter = std::make_shared<Dune::VTKWriter<GridView>>(gridView);
  VTKSequenceWriter<GridView> vtkSequenceWriter(vtkWriter, "getting-started-transport-fv-result");
  // Write the initial values
  vtkWriter->addCellData(c,"concentration");
  vtkSequenceWriter.write(0.0);
  // {construct_vtk_writer_end}

  // Inflow boundary values
  auto inflow = [](const FieldVector<double,dim>& x)
                {
                  return 0.0;
                };
  // Velocity filed
  auto v = [](const FieldVector<double,dim>& x)
           {
             return FieldVector<double,dim>(1.0);
           };

  // Now do the time steps
  // {time_loop_begin}
  double t=0;                  // Initial time
  const double tend=0.6;       // Final time
  const double dt=0.006;        // Time step size
  int k=0;                     // Time step counter

  while (t<tend)
    {
      // Apply finite volume schme
      paraevolve(gridView,mapper,c,dt,v,inflow);
      // Augment time and time step counter
      t += dt;
      ++k;
      // Write data. We do not have to call addCellData again!
      vtkSequenceWriter.write(t);
      // Print iteration number, time, and tim step size
      std::cout << "k=" << k << " t=" << t << std::endl;
    };
  // {time_loop_end}

  return 0;
}
